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In this report three types of probability distributions associated 
with random search-type models are derived. The distributions are closely 
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hunting or one-kind mine sweeping. The third type is multivariate allowing 
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1. Introduction 



In this report we investigate some discrete probability models which 
may be useful in search and related situations. The models are closely 
related to the classical occupancy problem of probability theory and may in 
fact be viewed as generalizations of it. 

The models are di screte in the sense that there are m discrete 
locations - referred to as cel 1 s - possibly containing the objects of the 
search - referred to as bal Is . The search itself is random in the sense 
that the cell to be searched is selected at random, i.e. with equal 
probability among all the cells involved. 

Two univariate models and a multivariate extension of the first of the 
two are considered. The detailed description of the first univariate 
model - called here Random Search-Variant 1 - is as follows: 

There are m cells, y Q of them initially empty and y^ of them 
containing a single ball each. The search consist of repeated independent 
trials, where, at each trial, a cell is selected with probability 1/m . The 
selected cell is examined, and if it contains a ball the ball is found with 
probability p . If a ball is found, it is removed from the cell so that 
the cell is now empty. Since the trials leading to a cell selection are 
assumed independent the same cell can be selected again regardless of the 
results of previous trials. Thus a cell can be selected even if a ball was 
found and removed from it at some previous trial. 

This model is investigated in Section 3. We find the probability dis- 
tribution of the number of balls found in n trials, the number of trials 
needed to find the k-th ball as well as moments of these distributions. 

As mentioned earlier, this model can be rephrased as an occupancy 
problem of distributing n balls into m cells. In fact, if y Q = 0 and 
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p = 1 it is identical with the classical occupancy problem - merely iden- 
tify the number of balls found in our model with the final number of nonempty 
cells in the occupancy problem. 

The second univariate model - Random Search Variant 2 - differs from 
Variant 1 only in the rule for selecting cells at each trial. This time 
we assume that the cell in which a ball has been found and removed is no 
longer a candidate for further selections. Thus, the number of cells to be 
searched is reduced by one each time a ball is found. The corresponding 
distributions of the number of balls found, the number of trials needed to 
find the k-th ball and their moments are described in Section 4. It may be 
interesting to notice that, with p = 1, this model is a hybrid of a classi- 
cal sampling with and without replacement. Think of the m cells as m 
balls, y Q of them white and y^ black. Then draw n balls replacing 
the white ones but not replacing the black ones, and ask about the 
distribution of the number of black balls drawn. 

The last section of this report deals with a multivariate extension 
of the Variant 1 model. The extension consists in assuming that there are 
now s different types of balls involved and each time a ball of type 
j > 0 is found it is replaced by a ball of type j - 1 or removed if 
j - 1 = 0 . The model parameters are now nonnegative integers y Q ,..., 
y s , where y Q +...+ y g = m , y Q is the number of empty cells and y. is 
the number of cells with balls of type j each initially. For simplicity, 
only the value p = 1 is considered here, i.e. if a cell containing a ball 
is selected the ball is found (and its type reduced) with probability one. 

The cell selection process is the same as in Variant 1. That is, cells 
are selected at random among the m cells in independent trials so that the 
same cell can be selected repeatedly regardless of the ball type it contains. 
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Again, this model can be rephrased as an occupancy problem of distributing 
n balls into m cells and interpreting a cell with ball type j as a cell 
containing s - j balls for j > 0 or a cell containing at least s balls 
for j = 0 . However, classical occupancy problems invariably assume that 
that initially all cells are empty, i.e. that y s = m in our notation. 
Therefore, classical results are not directly applicable to our model. 

In Section 5 we find the first two (joint) moments of the resulting 
numbers of balls type j = Q,...,s after n trials as well as of the num- 
bers of each ball type found. An expression for the marginal distribution 
of the number of balls type j is also found - although the complexity of 
this expression leaves some doubts about its possible uses. 

In this report, no attempt is made to study various asymptotic distri- 
bution resulting in letting the parameters of the models increase to infin- 
ity in some way. Although this may yield a considerable simplification, for 
instance by a Poisson distribution in the case of Variant 1, in most of the 
applications intended here the values of the parameters are typically small. 
Still, asymptotic results may perhaps warrant future investigation. 

Likewise, no algorithms to actually compute numerical values for the 
derived quantities are presented in this report - mainly because of a limited 
time and each of adequate computational facilities during preparation of 
this report. This, again, is left for a possible future work. 



3 



2. Possible Applications 



Although the purpose of this report is to develop model building tools 
rather then real world models themselves, a few possibilities will be 
mentioned briefly. 

A rather straightforward application of the two univariate cases would 
apply when one indeed has m distinct geographical locations out of which 
y-j contain at most one target each. Searching these locations at random 
would imply an absence of any particular search strategy or any further in- 
formation about the locations of these targets. Such a model may, for in- 
stance, serve as a basis against which an efficiency of various search 
strategies can be compared. In the case when the targets are distributed 
in a continuous geographical region, the discrete cells may be defined arti 
ficially by partitioning the region. The actual size of a cell must be 
chosen to satisfy the conflicting requirements of the size of the searcher' 
detection region and the assumption of a single target in each cell. This, 
of course, is the problem encountered any time a continuous model is being 
discretized. The choice between Variants 1 and 2 is then dictated by the 
nature of the targets and/or searcher. For example, if the targets are sta 
tionary. Variant 2 is more appropriate. Typical example here could be 
underwater mine hunting. If the targets are mobile and can randomly re- 
arrange between each search trial, Variant 1 could be used. Another possi- 
bility leading to Variant 1 is when the searcher (in the general sense of 
this term) has no control over the cell selection. Typical application in 
this case is of course an artillery coverage model - delivering n shells 
to a region partitioned into m cells out of which y - 1 contain a target. 
Similar applications, not involving an actual physical search, are con- 
ceivable. For instance, in reliability one can have an equipment with two 
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kinds of components - one with potentially unlimited supply of spares and 
the other with any a small number y-j of spares. 

The possible applications of the multivariate version are even more 
versatile depending on the interpretation of the ball type. Consider, for 
instance, a mine sweeping model with a ball type being the ship count the 
mine (ball) is initially set at. Alternately, consider a mine field with 
cells corresponding to paths through the mine field and ball type corre- 
sponding to the number of mines along a path. The cells would be separated 
from each other by the actuation width of a mine in this application. Yet 
another application is again the artillery coverage with multiple hits 
required to destroy a target. 

Our final remark concerns the choice of the model parameters. They 
can, of course, be left as true parameters to investigate their influence 
on the quantity of interest. But they themselves can be considered random 
variables resulting in various mixtures of the probability distributions 
involved. For instance, the number y-| of targets could be a Binomial ran- 
dom variable thus modeling the effect of a random appearance of targets. 

Even more realistically, perhaps, the ship counts in the mine sweeping ap- 
plication can be considered set at random among some given range of ship 
counts. Such assumptions could sometimes even simplify the analysis as 
shown e.g. at the end of Section 5. 



3. Random Search - Variant 1 



For the sake of reference let us restate the assumption of the model. 

We have m cells, y Q of them empty and y. j of them containing a single 
ball each. Cells are being selected repeatedly with equal probability 1/m 
in independent trials. If a selected cell contains a ball the ball is re- 
moved with probability p > 0 and remains there with probability q = 1 - p 
independently of the results of previous trials. Thus, in Variant 1, the 
same cell can be selected again regardless of the event that a ball may have 
been removed from it before. 

The parameters of the model are 

y 0 » y- 1 » and p , 

where y Q , y. | are nonnegative integers and 0 < p £ 1 . The letter 

m ■ y 0 * y i • 
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Let K( n ) be the number of balls obtained at the conclusion of the 
n-th trial. Clearly 0 £ K(n) <_min{n,y^} and 



(3.1) 



r 

K(n) 



K(n+1 ) = { 



K(n) + 1 



with probability 1 



y 1 - K(n) 



with probability 



y 1 - K(n) 



for n = 0,1,... with K(o) = 0 . 

Hence call ing 

P n (k) = P(K(n) = k) , 

we have the recurrence 



(3.2) 



p „ + i (k) 



m 



p n o<) 



^i 



- k + 1 



P P n (k-1) 



n = 0,1,... ; with P n (-1) = 0 and P Q (0) = ^ » P (k) = ^ ^ or ^ > 0 . 

This can be used to recursively evaluate the distributions p n (k) for 
given values of the parameters m, y-| and p . However, it is also 
possible to obtain an explicit formula, for instance by employing the 
generating functions 

oo 

(3.3) G. (t) = l t n P(k) , k = 0,1,. ..,y, . 

K n=0 n 1 
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From (3.2) with k = 0 we have 



W 0 * - 



yi 

1 -t" 



P (0) 
rr ' 



and using Pq(0) = 1 



P (0) = 
n v ' 



y l 

1 - — j- p 
m 



, n = 0,1 ,, 



Hence 



(3.4) 



G 0 (t) = 



1 



1 - t 




P 

; 



n+ i 

From (3.2) with k > 0 by multiplying by t n and summing over n = 0,1,. 
we obtain 



G k (t) 



1 - 



Yi - k 



G k (t) 



+ t 



yi - k + i 



m 



P G k _-|(t) 



whence solving for G k (t) and iterating 



(3.5) 



G k (t) k r 



y l 



(k) 



£t 

m 



n 

j=ol 



l - t 



y l ‘ J 

1 - -4— p 



m 



, k = 0,. . . ,y 1 



(k) y l J -1 

where yj = jy — . Since the k + 1 roots t^ = a^ , 

y l " J ' 

a. = 1 — p , of the denominator are all distinct we can use the 



partial fraction expansion 
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k 

n (l 

i=0 



ta.) 



I A.(l - ta.) 
j=0 J J 



-1 



where 



A. 

J 



a k 

J 



k 

n (a, 
i=0 J 



a i } 



-1 



j = 0,. . . , k 



y l ' J 

With a^ = 1 — p the coefficients are 



a j 



yi - j i 
i - 1 m p 

m r j 


fm) 

IpJ 


k 


k 


> 




n (j-i 
i=0 


) 













-1 



y i - J 'i 

1 n 


k 


V 


k 


m p j 






j i C k-j ) i 



Upon substitution into (3.5) we have 



G k (t) = yj k) t k I (-l) k_J 
K j=0 



y-, - j 

l - ~ p 



m 

j 1 ( k- j ) i 



(1 - ta.)' 1 , 



whence expanding (1 - ta.) - ^ into a geometric series 

sJ 



G k ( t) 



l t 

v=0 



v+k 



y l 

k 



l (-D 

j=0 



k-j fk 1 

J. 



1 - 



y! - j 



v+k 



m 
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The coefficients of this expansion are the desired probabilities, i.e. 



(3.6) 



P n (k) = 



y l 

i k ; 



l (-l) k ' j 
j=0 



w 



y-i - j ^ 



- \s. 

m 



n ryn k i 
1 11 i (-i ) j 

j-0 



k 

UU 



F " y l + k • j 



IP 



k , n — 0 , 1 ,... . 



Remark: 



(3.7) Since £ (-1 ) J ’ ( •] (x + k - j ) n = A k x n 

j=0 lJJ 



with a the forward difference operator one can also write 





fpl 


n 


r y O 


A k 


fm ) 




llj 




ik J 


A 


Ip • y 'J 



With y-| = m and p = 1 this is the classical occupancy distribution 
([FI], p. 58) although the generalizations with y^ < m or p < 1 can 
also be found in the literature ([JK], p. 124, 140). Note that the formula 
(3.6) gives automatically P (k) = 0 for k > n . 

Next, let us look at the moments of this distribution. The expectation 
u n = E[K(n)] (and higher moments as well) can be derived from the recurrence 
(3.2) by conditioning. For instance taking conditional expectation we 
obtain 
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n+1 



= (1 - £)u„ + 



nr "n 



y i 



£ 

m 



n = 0,1 ... . 



Uq = 0 , whence 



p n = y l 



1 - (i - £) 
m 



n = 0,1 , 



However, one can obtain all factorial moments 

= E[K(n)(K(n) - l)...(K(n) - r + 1)] , 

r = 1 . ,y-j , 
directly from (3.7) since 



y 



(r) = 
n 




v A k,m w \ n 

Cy-j-k) !(k-r) ! y l 

A j+r/m \n 
r) ^ ^1 ^ 



i ! 



by using the identity 



V 



l 

j=0 



V 




n 



(x * v) n 



n 



and linearity of the difference operator. Equivalently 



M 



(r) 

n 




I (-uMki -£j) 

j=0 J m 



r — 1 » • . • j y i > 

whence again the expectation 



‘n 1 ’ = ¥> - (1 -£>"l ' 



, ( 2 ) ( 2 ) 
From - y| ' 



1 - 2(1 - £)" H . (1 - 2 &)" 

m m 



we can get the variance 



* Var[K(n)J , 



^y,(l -£)"[i-y, O-f)"] 



. y p> tl . &>" . 



n = 0 ,1 , . . . , 



( 2 ) 

where y-j = y- ] (y 1 - 1) 
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Waiting time distributions: 



Consider next the random variable defined to be the trial number 

at which the k-th ball is found. We have 

W k T"i + . . . + , k - 1 , . . . ,y^ , 

where T. is the number of trials to find the j-th ball counted from the 

J 

the trial at which the (j-l)-th ball was found. Clearly, 

y-. - J + 1 

T. is geometric with the parameter - p and T-j , T^,... are 

independent. Thus calling 



>.(t) = l t n P(W. = n) 
K n=0 K 



the probability generating function of W^ we have immediately 



7] - J + 1 

k t — — = p 



m 



’k (t) ' " — y, - J + 1 

J i - tO - -J— p) 



m 



/ i \ a. k k-1 

y‘ k) 0 n 



j=0 



1 - t(l - 



7i - J 



vl 



m 



P) 



= (7i - k + D G k _-,(t) , 



where 6 is the generating function (3.3). 
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Hence 



^U) = l t v+k (y 1 - k + 1 ) £• 



;=0 



( v 'i u_i ( 

[y ’ ' V H) k -'- J ( k :') 

j=0 J 



1 

k-1 



1 - 



y l ~ J 



1 



v+k-1 



m 



and thus 



P(W k = n) = £ 



*1 

k 



j=0 1 J 



1 - 



*1 - J 



m 



n-1 



- fe] n k 



V k-i 



i (-d j 

j=0 



fk-1 

J J 



IF - », + k - 1 - J 



n-1 



n = k, k+1 ; k = 1 . ,y-| 



Again, this can be reexpressed as 



P(W k = n) = 



£l" k 



k-1 



J 



m 



*1 



n-1 



The expectation and the variance can be gotten most easily by using 
the fact that W. is a sum of independent geometric random variables T. 

K J 

wi th 



ElT.-i - Blv ”-+11 ■ Var[T.] = — * — j 

J P(y, J+D 0 p 2 (y r j+l) Z 



m 



y T J +1 

1 - -kr— P 



m 
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! 



so that 



and 



E ^] 



k-l 



p j=o y r J 



m k - 1 
Var[W ] • 2 i 
" P j=0 



(y-|-j) 



+ j 
2 
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4. Random Search-Variant 2 



This model differs from Variant 1 only in that the cell from which a 
ball has been removed is not selected again. The development parallels 
that of Variant 1 very closely and the same notation and definitions as in 
Section 3 are used. 

The basic recurrence for the random variables K(n) , n = 0,1,..., 
is now 



with K(0) = 0 . To avoid the trivial case yg = 0 we assume that y^ > 0 , 
i .e. y-| < m . 

The recurrence for the distribution P n (k) is now 



K(n) with probability 1 - 



Y) - K(n) 



m - K(n) P ’ 



K(n+1 ) = • 



K(n) + 1 with probability — v i \ p , 

m i\^ n j 



y-| - K(n) 




+ m - k + 1 



y-j - k + 1 



P P n (k-D , 



whence the generating functions G^(t) ar e 




G k (t) = IT 
n 

j=o 



16 



and calling 



y l " J 

a - = 1 - — — t- p 
J m - J v 



we have 



y Q p(j-i) 

3 j " a i (m-i ) (m-j ) * 



Thus 



A 



j 



y l - j 
m-j 




m - J 
(y 0 p) k 



k 

n 

i~0 

i7j 



m 

j 



2 

i 



1 



y l ~ J 
m-j 






P 

J 



m (k) (-D k - j 

(y 0 P) k 



and hence 



G k U) 



I ^ 4 

y 0 






v=0 



v+k 



l H> 

j=0 



k-j 



k 

Jj 



m - j 



From here the probabilities are 



p n^ "IT 



y l 

k 



y Ql 



' k 

l M) 

j=0 



k-j 



1 - 



y l - J 
m-j 



k 0 > • ♦ ♦ -j 5 n 0)1 9 ... • 
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Since with q = 1 - p 



1 



y-i - j 
m - j 



p = q + 



_ypP_ 

m - j 



we can also write 



p n ( k ) 



y P 



<V 1 
y l 



l (-D j 

j=0 



q + 



V 

m - k + j 






y 0 



r y ' 



& k F n (0) 



by using the identity 



wi th 



A k F(x) 



k I 

I (-1)' 

j=0 



F(x + k - j) 



F(x) 



q + 



V 

m - x 



The last expression can be used to obtain the factorial moments although 
the results are somewhat more complex that for Variant 1. We have for 
r - P,. . . ,y 1 



p 



(r) = 

n 



y i 

i 

k=r 



k (r) P n (k) 



y{^ y l- 

y£ +r j o 



V-j -r" 

. J . 



y^ A J (A r F n (P)) 
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by using the identity 



y(r) y r r 

n+r .^ n 
y 0 J=0 



l J ']' r ](y 0 - , ) J A rFn ( n -j) 



l fi)a b "^A^g(x) l "l(o-1) J g(x+b-j) 

j=0 j=0 



with b = y-j - r , a - and g(x) = A r F n (x) 



Removing the difference operator and substituting for F(x) we obtain the 
formula 



,(r) 



y, (r) r V , 

I l (->)' 



'1 
n+r 



q + 



y Q i=o j=o 

v | n 

m+i+j-n-r 





hi 


W 


. J . 



(yn-D j 



r = 0,...,y, 



In particular, the expectation is (r = 1) 



E[K(n)] = 




y r i 

l 

j=0 



fy r i) 



(yo - 1 ) 0 





f o + y ° P 1 


n 


r \ ~i 

„ + V ' 




q m-n+j-1 




^ m-n+jj 



The variance Var[K(n)] = + E[K(n)] - E“[K(n)] can be obtained 

similarly, resulting in a rather long expression. 
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Waiting time distribution - Variant 2 



The derivation is analogous to that of Variant 1 with the only difference 
•V] _ i + i 

being the parameter m ' j +" ' l P of the geometric random variables T. . 

The probability generating function ^(t) of the waiting time W^ 
is now 



so that 



y\ k) k k-1 
^(f) - / ry (pf) n 
K m [K> j=0 



1 - t 



1 - 



y l ~ J 
m - j 



y-i - k + i 

m - k + 1 pt G k-1^ 



V t v+k Pk 
v=0 y£(m-k+l) 



y l 

k 



k-1 

l (-D 

j-0 



k-1 -j fk-1 
. J 



- i 'l 



y l ~ J 
m - j 



v+k-1 



P(w k - n > - 



( \ 
y 



y Q (m-k+l) 



k-1 

l H) J 

j=0 



fk-1) 


/• ' 

O- 

o 

>> 

( 


> j J 


M m-k+l+j 



n = k, k+1 , . . . , k = 1 , . . . , y 1 . 
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From 



E 'V 



m - j + 1 

P(y r j+D ’ 



VartTj ] 



(m-j+ 1) 2 

p 2 (y-j - j + i ) 2 



y-, - j + i 

m - j + 1 



we get immediately 



and 



EI«kl * 7 l 



, k -1 m , 

I r rn - J 



P j = 0 y! - J 




k ; ] i_ 

p j=o *i - 



Var[W fc ] 



1 k ~ 1 

4 l 

P J =0 



•\ 2 r 



m-j 



y r j 




kq , y o k - v ] + ^0 

2 2 X n , ..2 

P p j=o (y-j -j ) 
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5. A Multivariate Search Model 



In this section we investigate some aspects of the multivariate gener- 
alization of the Variant 1 model. We have again m cells, yg of them 
empty, and m-yg of them containing one ball each. The balls can be of 
s _> 1 different types and initially, we have y-j cells containing balls 
of type 1, y 2 cells containing ball type 2, etc. We shall refer to the 
empty cells as containing ball type 0 , i.e. no ball. Thus the initial 
configuration is specified by a vector of nonnegative integers 

z - (y 0 » • • • » ) > 

where y~ + . . . + y = m and y. is the number of cells containing a ball 
US J 

of type j each. 

Cells are searched in repeated trials where at each trial a cell is 
selected at random equally likely among the m cells and independently 
of previous trials. If the selected cell contains a ball of type j > 0 
the ball is replaced by a ball of type j - 1 . Balls type 0 are just 
replaced, i.e. empty cells remain empty. The same cell can be selected 
repeatedly regardless of the result of previous trials. Thus for s = 1 
this is equivalent to the Variant 1 model with p = 1 . The model is 
therefore specified by a single vector parameter . 
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Let X_(n) = (Xg(n),..., X s (n)) T , n = 0,1,..., be random vector where 
for each j = 0, . . . , s 

x. t_ 

X.(n) = number of cells of type j after the n— trial. 

J 

Then _X(0) = and for n = 0,1,... 





X Q (n) 




with probability 


X-i (n) 

1 - . • 


X Q (n+l) = i 


X 0 (n) 


+ 1 


with probability 


X-j (n) 
m ’ 


(5.1) 


x j ( n ) 


- 1 


with probability 


X („) 
m ’ 


Xj (n+1 ) i 


Xj (n) 


+ 1 


with probability 


Xj+, ( n > 

m 




X . ( n ) 

J 




with probability 


Xj(n> ♦ X J+1 <n) 
m 


for 0 < j < s , and 










r 

X s (n) 


- 1 


with probability 


X s (n) 
m * 


X s (n+1) = i 


X s (n) 




with probability 


X (n) 

1 - -V- • 

m 



Of course 

(5.2) X Q (n) + X 1 (n) + ... + X $ (n) = m 
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for all n = 0,1 ,. . . . 
Alternately, let 



K(n) = (K Q (n),..., K $ (n)) 



be a random vector where for j = 0,..., s . 

K.(n) = number of trials (up to n-th) at which ball type j was found. 

J 

Then K(0) = 0 (zero vector) and 



K Q (n+l) = ^ 



Kg(n) + 1 with probability 
K Q (n) with probability 1 



y 0 + ^(n) 



m 



Y 0 + K-|(n) 



m 



K.(n) + 1 with probability 
J 



Kj(n+1 ) = • 



lV n) 



y, ♦ K i+1 (n) - K,(n) 



m 



with probability 1 



yj + K 1+ ,(n) - K.(n) 



m 



for 0 < j < s , and 



K (n) + 1 with probability 



K g (n+1) = - 



K s (n) 



with probability 1 



- K s (n) 
m 



y s ~ K s (n) 
m 



This time for all n = 0,1,... 



Kg(n) + ... + K s (n) = n 
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There is a simple linear relation between Xjn) and j<(n) , namely 



X Q (n) = y 0 + <i (n) , 



( 5 . 3 ) 



= y j ‘ l V n ) + K j+i^ n ^ ’ 0 < j < s » 



X s (n) - y s - K s (n) , 



and conversely for 0 < j <_ s 




wi th 



K Q (n) = n - (Kj (n) + . . . + K g ( n ) ) . 



Note that in terms of the mine sweeping model K^(n) is the number of 
deactivated (exploded) mines. Also B(n) = K-j(n) + ... + K g (n) , the 
total number of balls found, is the total number of contacts the sweeper 
makes during n sweeps. 

Either of the recurrence relations yields immediately the corre- 
sponding recurrence for the joint probability distribution of the random 
vector involved. For instance denoting 




• • • 9 



k $ ) = P(K 1 (n) = k-j , . . . ,K $ (n) = k $ ) 



we have for n = 0,1 



j • • • 
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J • • • ) 





S y' y .i + k .W - k 1 + 1 

£ m 



m 




m 




with initial condition Pg(k-| , . . . ,k s ) = 1 if k-j = . . . = k $ = 0 , and 
boundary conditions P n ( k-| , . . . ,k ) = 0 whenever k < 0 for some 
j = 1,.. . ,s . 

In principle, this can be used to recursively evaluate the distribution 
P (k-|,...,k ) but for larger values of s this is not practical since the 
storage requirements increase exponentially with s . Generally, an array 
of dimension (y $ + 1 ) (y $ + y $ _ -| + 1) ... (y + ... + y-| + 1) would be 
needed to store the current values of P n although some savings can still 
be made due to the fact that P„ = 0 unless 



We therefore attempt first to evaluate the first two moments of the 
random vector _X(n) . 

Define the multivariate moment generating function 



n 



°i k j + k j+i • 



s 



<» n (U = E exp l t £ X A (n) , 



£=0 



where t_ = (tg,...,t ) . From the recurrence (5.1) 
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s n 

exp l t X (n+1 ) = - 

_n £ £ m 



£=0 



X Q (n) 



+ l X (n)e 1-1 t? ' 
£=0 x 



exp l t X (n) . 
£=0 * 



Calling temporarily 



fl 

m 



»£<*> 



= ( 



ieVr^ 

m 



if 



if 



£ = 0 , 



0 < £ < s , 



and taking expectation 



( 5 . 4 ) 



Vl ^ 




9 

3t £ 



lp ( t ) 



Now 



9 Vi ( ^ 

at, 



s 

£=1 



3a, (t) 3^ n (t) 

3t i 3t £ 







a £ (t) 



3 2 ^ n (t) 

at .at ’ 

1 £ 



whence by setting t_ = 0_ . 
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S 3a £ (0) 



E[X,(n+l)] = l -rf E [X (n)] 

£=1 dt i 16 



+ E I EEX.(n)X (n)] 
m £-0 



s s 

However £ E[X • (n)X. (n) ] = E[X.(n) l X_(n)] = m E[X- (n) ] 
£=0 1 1 1 £=0 * 



since £ X (n) = m . 
£=0 * 



Further for 0 < £ < s 



f-J 

m 



3a £ (0) 

3t i"~ 



— if i = £ , 






if i = £ - 1 , 



otherwise , 



so that 



E[X 0 (n+l)] = E[X Q (n)] + i E£X ] (n) 3 



E[X.(n+l)J = (1 -l)E[X.(n)] + ^ E[X i+] (n) ] 



for 0 < i < s , and 



E[X s (n+l)] = (1 - J j r )E[X s (n)] 



Denoting X_(n) the column subvector 



(X-j (n) , . . . ,X s (n)) 
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and by Q the s * s matrix 



Q = 




] 



with 



we can write 



whence 



'u 



1 - - 
m 



*m 



if 



if 



i = J , 



j = i + 1 , 



otherwise , 



E[X(n+l)] = QE[X(n)] , 



E[X(n)] = Q n J 



with J = (y-, .... »y $ ) 



However, Q n is an upper triangular matrix with entries 



m 



n 



(m-1 ) 



n-j+i 



,(n) . 



= 1 






if 



if 



j > i , 



j < i , 



as is readily verified by induction. Hence 
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£[*,<"> 1 = -7t l j-i U"-i >"' J+1 y.i 

m j=i ^ J J 



(5.6) 



= \% (;](^) n ‘ £ y 

m £=0 W 



£+i ’ 



i = l,...,s , is the desired expectation vector. 



The expectation of X^Cn) is obtained from (5.2) 



E[X 0 (n)] = m - \ l l 
u m n i=l a=0 



S S-l 



n 

W 



Vi 



1 V 

= m - — l 

m n £=0 



s-l 



n 



(m-l )"- 1 Z, + , , 



where Z = y + . . . + y , £ = 1 ,. . . ,s . 

The expectation of the random vector K_( n ) is obtained using the linear 
relation (5.3), 



i s -' ' 



EIK (n)] - Z. l 

3 J m i=0 



"J Wi 



for 0 < j < s , and 



E[ Ko (n)] -Z, --L. I 
m £=0 



i s-l ( ^ 

1 r n 

M 



(m-1 ) n_£ Z 



£+1 • 
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Recurrence equations for the second moments can be developed similarly. 
Taking second partial derivatives of (5.4) we get 



s2 Vl(i) S 3 2 o,,(t) 3«. n (t) 



3t • 3t • 

* J 



S 3a ft) 3> n (t) 

+ y - 

3t i n j 3t t 



S Sa At) 3 * n (t) 

+ y — — 

A 3t • 3t-3t„ 

£=1 J It 



S 3 <|> (t) 

+ E 3t ~ at at 

z=0 z 3t i 3 j 3t £ 



Again upon setting _t = 0 the last terms becomes just E[X. (n)X.(n)] and 

J 

using (5.5) and 



3 \( 0 ) 

3t^ 3t . 



f 1 



' m 



if i = j = sl or i = j = z - 1 , 



otherwise , 



where 0 < i < j £ s 



0 < i < s , we obtain for 0 < i < j <_ s 
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E [ X _j (n + 1) Xj (n + 1) ] = E[X i (n )X - (n ) ] 

^E[X f (n)X 1+ i(n)] + i E[ x .(„)x j+ i(n)] 

! S ECX i (n)] + H E[x m (n)] if J- '' 

- ^ E [ X j ( n ) ] if j = i + 1 , 

0 if j > i + 1 , 

with terms E[X s (n )X s+ ^(n )] and E[X $+ ^(n)] being interpreted as zero. 
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This again looks better in matrix notation. Denote 



M „ ■ • A - • E „ * 
the s x s matrices defined by 



u j ( n ) = E[X 1 (n)Xj(n)] , 



if 



J = i . 



a ■ ■ 
TJ 



if 



j = i + 1 , 



otherwise , 



m 



1 E[X ( (n)] +iE[X 1 + 1 (n)] 



E ij (n) 



iE[X j( n)] 

iE[X t (n)] 



if 

if 



J = 



if 

i + 1 , 
j + 1 , 



0 



otherwise. 



i , 



Then for n = 0,1,... 

(5.7) M n ., = M + AM + M A T + £ n , 

n+1 n n n n 

with T denoting transposition. 

Since the entries of Mq are just y^y^ while the entries of E n 
are available from (5.6) for all n = 0,1,..., the second moment matrix 
M can be evaluated iteratively. 
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It is even possible to write a closed form expression since upon iterating 
(5.7) we obtain 



(5.8) 



n v 

- I l 

v=0 1=0 



(")(>V A 



v-£ jT 



n-1 k v 

+ l l l 

k=0 v=0 1=0 






'n-l-k 



(A U ~V 



n = 1,2,... . The £-th power A = [a: . ] of the matrix A is an 

^ J 

upper triangular matrix with entries 




if j > i , 



if J < i , 



so that all the terms on the right-hand side of (5.8) are available. 

Having the mean vector and the second moment matrix one can calculate 
the covariance matrix, from which the covariance matrix of the random 
vector l<(n) is obtain by using the linear transformation (5.3). Although 
having an expression for the mean vector and the covariance matrix is use- 
ful one would like to obtain a formula for the probability distribution as 
we did in the univariate case. Unfortunately, the multivariate case is 
considerably more complicated and only a partial result is obtained. 
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Let 0 < j <_ s be fixed and denote by x(^(n) the number of cells 
among the original cells (with ball type i initially) that after n 

trials contain balls of type j . Clearly 



where of course 



Xj (n ) = Xj 0) (n ) + ... + X^. s) (n) , 

X^ (n ) = 0 for i < j . 

J 



Let us now condition on the random vector N_ = (N ,...N ) , where N . 

is the number of trials which resulted is a selection of a cell from the 

original group of y. cells. Then N_ is multinomial with parameters n 

and y./m , i = 0,...,s , and X^(n) ,..., x( s ^(n) are conditionally 

J j 

independent given N. . It follows that 



( 5 . 9 ) 



P(X J (n) = k) 





s 



m 



n 



P(x(°^(n) = '< | N_) * P(x(^(n) = k|N) *. . . 

J J 

* P(X< s) (n) = k|N) , 



where the summation is over all nonnegative integers n Q ,...,n s such that 

n + ... + n = n and asterisks denote convolutions, 
o s 

Now for i < j we have trivially 

P(xj lJ (n) = 0|N) = 1 , 
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while for j _< i £ s the conditional distri but ion P(x( 1 ^ (n ) = k |_N_ = n_) 
is identical with the classical occupancy distribution P^(k) of finding 
exactly k cell occupied with exactly i - j balls after distributing n^ 
balls randomly into y = y. initially empty cells. 

This distribution is most easily derived from the generating function 



H ( y ^ ( u , v ) 



- I l 

n=0 k=0 



(yu) 



v k P 



(y) 



OO 



by conditioning on the number of balls placed in each cell. (See [JK], 
p. 116 for details.) It follows that 

H^(u,v) = (u,v)j y 

and since P^(l) = 1 if n = i - j while P^(0) = 1 otherwise we have 
n n 

(5.10) H (y) (u,v) = 

J 



(v-1) 



u 1_J 

TmT T 



+ e 
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Now for 0 < k < y 



S k H</> 

J 


- Ml u f ‘ J \ 


u u'-J 


> 

CO 


■' 5 

1 

1 

O 

II 

> 


thtt 



y-k 






i-j k+£ 






(y-k-A)u 



(k) 



l "l t-D l ( y ; k ) (y-^) u - v 

= n o -n ' x, / *■ 



v=0 £=0 



v!((i-j)!) 



so that P^(k) = ^ x coefficient at u 

y 



n 



Thus with y = y. and n = n. we have 

J J l l 



(y i ) 

p 1 (k) 

n i n i t=0 



y- k) n i ! y i _k i -k \ 

1 1 i (-in \ ) 



y 



n i - (i-j)(k+£) 

(y-k+t) 

((1-j)!) k+ *(n. - (i - j ) (k+2.) ) ! * 



Upon substitution into (5.9) we obtain the desired expression for tne 
probabilities (for 0 < j <_ s ) 
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• « t 



P(X f (n) *k) -4 I 1 ^ 2 - 
J m 






<„j + 



+ n s> 



(n - (n. + ... + n s ))! 



y i' k i 

I 

£.=0 

1 



(-1)^ 

Z. ! (y.-k . ) ! 

1 1 l 1 



(y fW 



n 1 -(i -j ) (k 1 +a i ) 



k +£ . 

((i-j)!) 1 i (n.-(i-j)(k i +&.))! 



where the first two sums are, respectively, over all nonnegative integers 
n.,...,n and k.,...,k such that 

J ^ J ^ 

n . + ...+ n _< n , 

J ^ 

and 

k . + ... + k = k . 

J ^ 



Unfortunately, for larger values of s, k, or n the number of terms 
involved in the indicated sums is rather prohibitive. The problem of 
developing a manageable computational algorithm to evaluate the probabili- 
ties P(X.(n) = k) is left for future investigation. 

J 

It should be pointed out, however, that the generating function 

G,(u,v) = l l v k P( X . (n ) = k) 

J n=0 k=0 J 



is just a product of the generating functions (5.10), specifically 



Gj(u,v) = e 



(y n + -*- + y^ (y,-) 



(y e ) 



^ _1 Hj J " ( u , v ) . . . H j S (u,v) . 
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Thus, if an algorithm can be found to expand the right-hand side of this 
expression into a power series in u and v the probabilities could be 
determined directly from this expansion. 

As an example of this approach consider the case when the parameters 
y lt ...,y are themselves random, having a multinomial distribution with 
parameters s , b =yj+...+y s * and Pj = 1/s , i = l,...,s , respectively. 
In other words assume that each ball type 0 < j <_ s was selected independ- 
ently and at random among all nonzero ball types l,...,s . Their using the 
multinomial theorem we see that 



S 1 (u , v ) = e 



y u 
J 0 



s i-lf 

be “ + (v ' u 4i iTTirr 



Then taking k-th partial derivative with respect to v and setting 
v = 0 we get 



b 



(k) 



b-k 



l 

1=0 



(-D 



% 



s --> uM b (b-k-t) 
i =0 1 y 



(m-k-t)u 

e v 



After expanding the exponential function it is seen that the coeffi- 
cients of powers of u will be expressible as ordinary double sums, 
involving coefficients of powers of the polynomial 



The latter could possibly be precomputed and thus there is a hope for 
a reasonably efficient algorithm. 
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